FROZEN GAUSSIAN APPROXIMATION FOR HIGH 
FREQUENCY WAVE PROPAGATION 



JIANFENG LU AND XU YANG 

Abstract. We propose the frozen Gaussian approximation for computation 
of high frequency wave propagation. This method approximates the solution 
to the wave equation by an integral representation. It provides a highly ef- 
ficient computational tool based on the asymptotic analysis on phase plane. 
Compared to geometric optics, it provides a valid solution around caustics. 
Compared to the Gaussian beam method, it overcomes the drawback of beam 
spreading. We give several numerical examples to verify that the frozen Gauss- 
ian approximation performs well in the presence of caustics and when the 
Gaussian beam spreads. Moreover, it is observed numerically that the frozen 
Gaussian approximation exhibits better accuracy than the Gaussian beam 
method. 



1. Introduction 

Wc arc interested in developing efficient numerical methods for high frequency 
wave propagation. For simplicity and clarity we take the following linear scalar 
wave equation to present the idea, 

(1.1) d?u-c 2 {x)Au = 0, xeM. d , 

with WKB initial conditions, 



(1.2) 



u (x) = A Q (x)ei So< - x \ 
d t u (x) = ±B (x)e7 s ^\ 



where u is the wave field, d is the dimensionality and i = y— 1 is the imaginary 
unit. We assume that the local wave speed c(x) is a smooth function. The small 
parameter e <^ 1 characterizes the high frequency nature of the wave. The proposed 
method can be generalized to other types of wave equations |17j . 
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Numerical computation of high frequency wave propagation is an important 
problem arising in many applications, such as electromagnetic radiation and scat- 
tering, seismic and acoustic waves traveling, just to name a few. It is a two-scale 
problem. The large length scale comes from the characteristic size of the medium, 
while the small length scale is the wavelength. The disparity between the two length 
scales makes direct numerical computations extremely hard. In order to achieve 
accurate results, the mesh size has to be chosen comparable to the wavelength or 
even smaller. On the other hand, the domain size is large so that a huge number 
of grid points are needed. 

In order to compute efficiently high frequency wave propagation, algorithms 
based on asymptotic analysis have been developed. One of the most famous exam- 
ples is geometric optics. In the method, it is assumed that the solution has a form 
of 

(1.3) u{t,x) = A{t,x)e lS ^ x)/e . 

To the leading order, the phase function S(t, x) satisfies the eikonal equation, 

(1.4) \d t S\ 2 -c 2 (x)\V x S\ 2 =0, 
and the amplitude A(t, x) satisfies the transport equation, 

9, ,V X S (d 2 S -c 2 (x)AS) „ 

The merit of geometric optics is that it only solves the macroscopic quantities 
S(t, x) and A(t, x) which are e-independent. Computational methods based on the 
geometric optics are reviewed in [31126]. 

However, since the eikonal equation (|1.4p is of Hamilton- Jacobi type, the solution 
of (|1.4[) becomes singular after the formation of caustics. At caustics, the approx- 
imate solution of geometric optics is invalid since the amplitude A(t, x) blows up. 
To overcome this problem, Popov introduced Gaussian beam method in [19 . The 
single beam solution of the Gaussian beam method has a similar form to geometric 
optics, 

u(t,x) = A(t,y)e lS{t ' x ' v)/e . 
The difference lies in that the Gaussian beam method uses a complex phase function, 

(1.5) S(t, x, y) = S(t, y) + p(t, y)-(x — y) + -(x — y) ■ M (t, y)(x - y), 

where S G M, p G M. d , M £ C dxd . The imaginary part of M is chosen to be 
positive definite so that the solution decays exponentially away from x = y, where 
y is called the beam center. This makes the solution a Gaussian function, and 
hence the method was named the Gaussian beam method. If the initial wave is not 
in a form of single beam, one can approximate it by using a number of Gaussian 
beams. The validity of this construction at caustics was analyzed by Ralston in 
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[24] . Recently, there have been a series of numerical studies including both the 
Lagrangian type [Hl[22l[23j[28l[29] and the Eulerian type [EHUMHUHllinilS] - 

The construction of Gaussian beam approximation is based on the truncation 
of the Taylor expansion of S around the beam center y up to the quadratic term, 
hence it loses accuracy when the width of the beam becomes large, i.e., when the 
imaginary part of M(t, y) in (|1.5[) becomes small so that the Gaussian function 
is not localized any more. This happens for example when the solution of the 
wave equation spreads (the opposite situation of forming caustics). This is a severe 
problem in general, as shown by examples in Section [4] One could overcome the 
problem of spreading of beams by doing reinitialization once in a while, see |22p23j . 
This increases the computational complexity especially when beams spread quickly. 

Therefore a method working in both scenario of spreading and caustics is re- 
quired. The main idea of the method proposed in the current work is to use 
Gaussian functions with fixed widths, instead of using those that might spread 
over time, to approximate the wave solution. That is why this type of method is 
called frozen Gaussian approximation (FGA). Despite its superficial similarity with 
the Gaussian beam method (GBM), it is different at a fundamental level. FGA is 
based on phase plane analysis, while GBM is based on the asymptotic solution to a 
wave equation with Gaussian initial data. In FGA, the solution to the wave equa- 
tion is approximated by a superposition of Gaussian functions living in the phase 
space, and each function is not necessarily an asymptotic solution, while GBM uses 
Gaussian functions (named as beams) in the physical space, with each individual 
beam being an asymptotic solution to the wave equation. The main advantage of 
FGA over GBM is that the problem of beam spreading no longer existsQ Besides, 
numerically we observe that FGA has better accuracy than GBM when keeping the 
same order of terms in asymptotic series. On the other hand, the solution given 
by FGA is asymptotically accurate around caustics where geometric optics breaks 
down. 

Our work is motivated by the chemistry literature on the propagation of time 
dependent Schrodinger equation, where the spreading of solution is a common phe- 
nomenon, for example, in the dynamics of a free electron. In [5], Heller introduced 
frozen Gaussian wavepackets to deal with this issue, but it only worked for a short 
time propagation of order 0(H) where H is the Planck constant. To make it valid for 
longer time of order 0(1), Herman and Kluk proposed in [B] to change the weight 
of Gaussian packets by adding so-called Herman-Kluk prefactor. Integral represen- 
tation and higher order approximations were developed by Kay in [12] and [13]. 
Recently, the semiclassical approximation underlying the method was analyzed rig- 
orously by Swart and Rousse in [27] and also Robert in [25]. We generalize their 



Divergence is still an issue for the Lagrangian approach, one needs to work in the Eulerian 
framework to completely solve the problem, which is considered in |17| . 



4 



JIANFENG LU AND XU YANG 



ideas for propagation of high frequency waves, aiming at developing an efficient 
computational method. We decompose waves into several branches of propagation, 
and each of them is approximated using Gaussian functions on phase plane. Their 
centers follow different Hamiltonian dynamics for different branches. Their weight 
functions, which are analogous to the Herman-Kluk prefactor, satisfy new evolution 
equations derived from asymptotic analysis. 

The rest of paper is organized as follows. In Section^ we state the formulations 
and numerical algorithm of the frozen Gaussian approximation. In Section [3l we 
provide asymptotic analysis to justify the formulations introduced in Section[2j The 
numerical examples are given in Section|4]to verify the accuracy and to compare the 
frozen Gaussian approximation (FGA) with the Gaussian beam method (GBM). 
In Section [SJ we discuss the efficiency of FGA in comparison with GBM and higher 
order GBM, with some comments on the phenomenon of error cancellation, and we 
give some conclusive remarks in the end. 

2. Formulation and algorithm 

In this section we present the basic formulation and the main algorithm of the 
frozen Gaussian approximation (FGA), and leave the derivation to the next section. 

2.1. Formulation. FGA approximates the solution to the wave equation by 
the integral representation, 

u FGA (t,x) = 1 / a + (t, q ,p)e^ tx 'y^u + ,o(y)dydpdq 
, 2 ^ (2ire) 3d / 2 J R 3 d 



^ / a_(t,q,p)ei*-^y^u^.o(y)dydpdq, 



(27T£ 

where u±,o are determined by the initial value, 

(2.2) u ±i0 (x) = A ± (x)ei s °W, 

with 

A ± (x) = -[Ao(x)± 



2 V c(x)\d x S (x) 
The equation (|2.1[) implies that the solution consists of two branches ("±"). 
In (|2.ip . $± are the phase functions given by 

(2.3) $±(t,x,y,q,p) = P±(t,q,p) ■ (x - Q±(t,q,p)) - p ■ (y - q) 

+ \\x~Q ± {t,q,p)\ 2 + l -\y-q\ 2 . 

Given q and p as parameters, the evolution of Q ± and P± are given by the equation 
of motion corresponding to the Hamiltonian H± = ±c(Q±)\P±\, 



(2.4) 



^ = -d Q± H ± = T d Q± c\P ± \ 
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with the initial conditions Q±(0,q,p) — q and P±(0,q,p) = p. The evolution 
equation of a± is given by 



^r = ± ^\p7\- dQ - c -lp7T c 

(2.5) ± ^J Z ^d z Q ± (2^-®d Q± c 



'P ± ®P ± 



l)-^P ± \d 2 Q± cfj 



\P±\\ \P±\ 2 
with the initial condition, 

a ± (0,q,p) = 2 d / 2 . 

In (12.51) . P± and Q± are evaluated at (t, q,p), c and Oq ± c are evaluated at Q ± , / 
is the identity matrix, and we have introduced short hand notations 

(2.6) d z = d q - id p , Z ± = d z (Q ± + tP±). 

The evolution of the weight a± is analogous to the Herman-Kluk prefactor [6]. 

Remark. 1. The equation (|2.5|) can be reformulated as 

/„ - da± P± a± / 1 dZ±\ 

(2J) -dT = ±a± |Pi|-^ ± c+ T tr ^ —J- 

When c is constant, (|2.7[) has an analytical solution a± = (detZ-t) 1 / 2 with the 
branch of square root determined continuously in time by the initial value. 
2. d z Q± and d z P± satisfy the following evolution equations 

d(d z Q ± ) d Q± c®P ± 
( 2 - 8 ) T7 =±d z Q ± r— -j ±cd z P± 



dt ^ ± \P ± \ - z =\\P ± \ |P±| 3 

One can solve (|2 . 8|) - (|2.9() to get d z Q ± and d z P± in (|2.5[) . This increases the 
computational cost, but avoids the errors of using divided difference to approximate 
derivative. 



Notice that (|2.1j) can be rewritten as 

a+ 

t (27re) 3d / 2 



(2.10) 



where 



u FGA (t,x)= / - a + /9 ^ +e ^ p +-^-Q + )-Al^-Q + l 2 d pdg 



^_ e JP_.( a ._Q_)_i|x-Q_|- dpdg ^ 



>. d (2-rre) 3d /< 



(2.11) i>±(q,p)= [ u ± , (y)e-iP<y-^-^y-rf dy. 

Therefore, the method first decomposes the initial wave into several Gaussian func- 
tions in phase space, and then propagate the center of each function along the 
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characteristic lines while keeping the width of the Gaussian fixed. This vividly 
explains the name frozen Gaussian approximation of this method. 

The formulation above gives the leading order frozen Gaussian approximation 
with an error of 0(e). It is not hard to obtain higher order approximations by 
the asymptotics presented in Section [3] We will focus mainly on the leading order 
approximation in this paper and leave the higher order corrections and rigorous 
numerical analysis to future works. 

2.2. Algorithm. We first give a description of the overall algorithm. To construct 
the frozen Gaussian approximation on a mesh of x, one needs to compute the 
integral (|2 . 10[) numerically with a mesh of (q,p). This will relate to the numerical 
computation of p. lip with a mesh of y. Hence three different meshes are needed 
in the algorithm. Moreover, the stationary phase approximation implies that tp± in 
(|2.1ip is localized around the submanifold p — V q 5o(q) on phase plane for WKB 
initial conditions (II. 21) when e is small. This means we only need to put the mesh 
grids of p around V q So(q) initially to get a good approximation of the initial value. 
A one-dimensional example is given to illustrate this localization property of ip± in 
Figured] (left). The associated mesh grids are shown in Figured] (right). 




q q 



Figure 1. Left: an illustration of the localization of tp+ on (q,p) 
domain for u +i o(y) = exp ^z sl "^^ , e = 1/128; the black solid 
curve is p = cos(6g)/2. Right: the corresponding mesh grids of 

Next we describe in details all the meshes used in the algorithm. 

(1) Discrete mesh of (q,p) for initializing Q,P. Denote 8q — (Sqi, ■ ■ ■ ,Sqd) 
and Sp = (Spi, • • • , 5pd) as the mesh size. Suppose q° = (g°, • • • , gj]) is the 
starting point, then the mesh grids q k , k — (fci, • • • , fed), are defined as 

q k = {q°i + (fci - !)£<?!, ' • • ,«S + (kd ~ l)Sq d ), 
where kj = 1, • • • , N q for each j £ {!,■■■ , d}. 
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The mesh grids p k ' e , £ = (ii, ■ ■ ■ ,£d), ar e defined associated with the 
mesh grids q k , 

p k '*=(d qi S (q k ) + h6 Pl ,--- ,d qd s (q k )+e d 5p d ), 

where ij — —N p , ■ ■ ■ ,N P for each j G {1, • • ■ , d}. 

(2) Discrete mesh of y for evaluating ip± in (|2.11j) . 8y = {6y\, • • ■ , 5yd) is the 
mesh size. Denote y° — (yj, • • • ,y d ) a s the starting point. The mesh grids 
y m are, m = (mi, • • • ,m d ), 

V m = (?/? + (mi - 1)<W • • ,y° d + {md-l)5yd), 

where rrij = 1 , • • • , N y for each j £ {1 , • • • , d} . 

(3) Discrete mesh of x for reconstructing the final solution. Sx = {8x\ , • • • , 8x d ) 
is the mesh size. Denote a; = [x\, • • • , x d ) as the starting point. The mesh 
grids x n are, n = (m, • • • , n d ), 

x n = (xl + (ni - l)6x u • ■■ ,x° d + (n d - l)Sx d ), 

where rij = 1, ■ ■ ■ , N x for each j £ { 1 , ■ • • , d} . 

With the preparation of the meshes, we introduce the algorithm as follows. 

Step 1. Decompose the initial conditions (|1.2[) into two branches of waves according 
to (I2T21) . 

Step 2. Compute the weight function ij)± by (|2.1ip for (Q, P) initialized at (q k ,p k ' £ ), 

(2.12) ^ ± {q k ,p k ' 1 ) = ^eK-^'^-J^+fl^-^l 2 ) 
m 

x u ±fi {y m )re{\y m - q k \)5yi ■ ■ ■ 6y d , 

where re is a cutoff function such that rg = 1 in the ball of radius 9 > 
centered at origin and re = outside the ball. 
Step 3. Solve (|2.4I) - (|2.5[) with the initial conditions 

Q ± (0,q\p k > e ) = q k , P ± (0,q k ,p k < e )=p k >\ 
a±(0,q fc ,p fc ' £ )-2 d / 2 , 

by standard numerical integrator for ODE, for example the fourth-order 

Runge-Kutta scheme. Denote the numerical solutions as (Q± £ , P± ) and 
k,e 
a ± ■ 

Step 4. Reconstruct the solution by (|2.10[) . 



FGA (t,x n ) = y( a+ 2M l k .p k *)e^ p Y-^-QY)-h\^-QY\ 

(2> ^_ iq k ip k,ty iP ^-( x "-Q^)-±\ x ™-Q h jr 



(27T£) M / 2 ' 

x 5gi • • • <5q d (5pi • • • 6 P d, 
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where r g = rg(\x 




Remark. 1. In setting up the meshes, we assume that the initial condition (|1.2j) 
either has compact support or decays sufficiently fast to zero as x — > oo so that we 
only need finite number of mesh points in physical space. 

2. The role of the truncation function rg is to save computational cost, since 
although a Gaussian function is not localized, it decays quickly away from the 
center. In practice we take 9 = 0(- s /e), the same order as the width of each 
Gaussian, when we evaluate (|2.12l) and (|2.13p numerically. 

3. There are two types of errors present in the method. The first type comes from 
the asymptotic approximation to the wave equation. This error can not be reduced 
unless one includes higher order corrections. The other type is the numerical error 
which comes from two sources: one is from the ODE numerical integrator; the 
other is from the discrete approximations of integrals (|2.10l) and (|2.1ip . It can 
be reduced by either taking small mesh size and time step or using higher order 
numerical methods. 

4. Note that the assumption that the initial conditions are either compactly 
supported or decay quickly implies that the values on the boundary are zero (or 
close to zero). Then (|2.12j) and (|2.13p are the trapezoidal rules to approximate 
(12.111) and (|2.10l) . Notice that, due to the Gaussian factor, the integrand functions 
in (|2.11[) and (|2.10[) are exponentially small unless x — Q and y — q are in the 
order of ©(e 1 / 2 ), which implies their derivatives with respect to y, q, p are of the 
order C(e -1 / 2 ). This suggests Sy, 6q, Sp should be taken as the size of 0(^/e). 
Hence N y and N q are of order 0(e~ d l 2 ). As illustrated in Figure [I] N p is usually 



taken as O ^ min ,; p J , which is of order 0(1). N x is not constrained by e, and is 
only determined by how well represented one wants the final solution. 

5. Step 2 and 4 can be expedited by making use of the discrete fast Gaussian 
transform, as in [2"2"ll23| . 

3. Asymptotic derivation 

We now derive the formulation shown in Section [5] using asymptotic analysis. 
We start with the following ansatz for the wave equation (jl.ll) . 



(3.2) $± (t, x,y,q,p) = S± (t, q, p) + P± (t, q, p) ■ {x - Q ± (t, q,p)) - p ■ {y - q) 



(3.1) 




w+,o(y)dydpd<7 




where $± are given by 



+ ^\x-Q ± (t,q,p)\ 2 + -\y-q\ 2 . 



FROZEN GAUSSIAN APPROXIMATION 9 

The initial conditions are taken as 

Q±{Q,q,p) = q, P±(o,q,p)=p, 

(3 3) 

S ± (0,q,p) = 0, a ± (0,q,p) = 2 d / 2 . 

The subscript ± indicates the two branches that correspond to two different 
Hamiltonian, 

(3.4) H + (Q + ,P+) = c(Q+)\P + \, H_(Q_,P_) = -c(Q_)|P_|. 
J 3 -!- and Q ± satisfy the equation of motion given by the Hamiltonian H± 

P± 



(3.5) 



d t Q± = d P± H± = ±c- 

fyP ± = -d Q± H ± = T d Q± c\P±\. 

By plugging (|3.1[) into (11.11) . the leading order terms show that the evolution of 
S± simply satisfies 

(3.6) d t S± = 0. 

This implies S±(t,q,p) = 0. Hence we omit the terms S± in Section [5] and later 
calculations. 

Before proceeding further, let us state some lemmas that will be used. 
Lemma 3.1. For u G L 2 (W l ), it holds 

(3.7) u(x) = 1 / 2 d /V*±(°™^u(y)dydpd<z. 



Proof. By the initial conditions (I3.3[) . 



l - '2 , „,2 



(3.8) $±(0, x,y,q,p) = p ■ (x - q) - p ■ (y - q) + -\x - q\ 2 + -\y - q 

Therefore, (|3 . T[) is just the standard wave packet decomposition in disguise (see for 
example [4]). □ 

The proof of the following important lemma follows the one of Lemma 3 in [27j . 

Lemma 3.2. For any vector a{y,q,p) and matrix M(y,q,p) in Schwartz class 
viewed as functions of (y,q,p), we have 

(3.9) a(y, q,p) ■ (x - Q) ed Zh (ajZj£), 

and 

(3.10) (x-Q)- M(y, q,p){x - Q) ~ ed^Q.M^ 1 + e 2 d Zm (d Zl (M jk Zrf)Z£), 

where Einstein's summation convention has been used. 
Moreover, for multi-index a that \a\ > 3, 

(3.11) (x-Q) a -©(e 1 " 1 " 1 ). 
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Here we use the notation f ~ g to mean that 



(3.12) 



/ /e^ $± dy dpdq = / ge?' s ' ± dy dpdq. 



Proof. Since the proof is exactly the same for the cases of $+ and we omit the 
subscript ± for simplicity. As a and M are in Schwartz class, all the manipulations 
below are justified. 

Observe that at t = 0, 

-{d q Q)P + p = 0, {d p Q)P = 0. 

Using (|3.5p . we have 

d t (-(d q Q)P + P ) = -d q (d t Q) P - d q Qd t P 

= -d q {c^jP + d q Qd Q c\P\ 
= 0. 

Analogously we have dt((d p Q)P) — 0. Therefore for all t > 0, 
-{d q Q)P + p = 0, {d p Q)P = Q. 
Then straightforward calculations yield 

9 q $ = (d q P - id q Q)(x Q) i(y - q), 
= (d p P - id p Q){x -Q)-(y- q), 

which implies that 

(3.13) id z <$> = Z(x - Q), 

where d z and Z are defined in (|2.6[) . Note that, Z can be rewritten as 

'd q Q d q P\ f-il\ 

q d p p) \ i J ' 

where / stands for the d x d identity matrix. Therefore, define 



Z = 



W + .P)=(,7 



F 



d q Q d q P\ 
d p Q dpP ) ' 



then 



■I I) FF T I ~ l J + 2F 
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In the last equality, we have used the fact that 

f(° -*V=(° 

\tl J \il o J 1 

due to the Hamiltonian flow structure. Therefore ZZ* is positive definite for all t, 
which implies Z is invertible and 

(3.14) (x-Q) =iZ~ 1 d z <f>. 

Using (|3.14p . one has 

/ a- (x-Q)ei*dydpdq = e [ a 3 Z^ (-d Xk $) dy dpdq 

jR3d J R 3d \e / 

= -e / d Zk (a j Z~ k 1 )e-' s ' dydpdq, 

where the last equality is obtained from integration by parts. This proves (|3.9j) . 
Making use of (|3.9j) twice produces (|3.10p 

(x-Q)- M(x -Q) = (x- Q)jM jk (x - Q) k 

^-ed Zl ((x-Q)jM jk Z^) 
= ed zl QjM jk Z^ - e(x - Q) 3 d Zl (M jk Z^) 
~ sd zl QjM jk Z^ + e 2 d Zm (d Zl (M jk Z^)Z^). 
By induction it is easy to see that (|3.11[) is true. 



□ 



3.1. Initial value decomposition. By (|3.2j) and (|3.5p we obtain that 

ft*± = P± ■ d t Q ± + {d t P ± - id t Q ± ) ■ (x - Q ± ) 

(3.15) t p, v 

= T c\P ± \T(x-Q ± )-(\P±\d Q± c + i-±-c), 



and in particular for t = 0, 

(3.16) d t $±(0, x,y,q,p)=Tc\p\T(x-q)- (\p\d q c + z^c) . 

The ansatz (|3.ip shows that 

u(0,x)= fo 1 / a+ (0,g,p)ei*+^*^ U+ . o (y)d 2/ dpd 9 



(3.17) 

and 
(3.18) 

3*u(0, sc) 



™ / a_(Q 1 p J g)e**-(°'»'»'«*)«_ > o(»)d»dpdii 1 



(27TE 



(0 \ 3d/2 [ {d t a + + l ^d t $ + )e^ <*^u + , (y)dydpdq 
' / (dta- + — a t $_)e^-^^ p ) U _.o(y)dydpdq. 



(2^e)3rf/2 y R3d 
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We take 

(3.19) a±{0,q,p) = 2 d '\ 

(3.20) «+,o(as) = A + (x)ei So(x \ 

(3.21) u-,o(as) = A_{x)e^ Sa{ - x \ 
with 

(3.22) A ± (x) = Ua (x)± ' LUx] 



2V ' ' c(x)|9 x ,So(x)|, 
We next show that this will approximate the initial condition to the leading order 
in £. 

Substituting (|3.19[) - (|3.22[) into (j3. 1 7[) and using Lemma I3.1[ we easily confirm 
that 

u(0, x) = u +fi {x) + u- fl (x) = A (x)ei So{x) . 

For the initial velocity, we substitute p.l9[) ~ p.22[) into (j3. 18[) and keep only the 
leading order terms in e. According to Lemma I3.2| only the term =pc|f>| in dt<&± 
will contribute to the leading order, since the other terms that contain (x — q) are 
0(e). Hence, 

2 d/2 



nd/2 r 

d t u(0,x)=- - / - C (g)|p|e**+<°'-*«*>u + ,o(y)dydpd* 

! -c(q)\ P \ei*-^y^U-, (y)dydpdq + O(l). 

jR3d e 



2 d/ 



(2ne) 3d / 2 
Consider the integral 

/ c(q)\p\e-^ p - {y - q) -^ y - q \ 2 A±{y)e L c So ^ &y = [ c(q)\p\A± (y)ei @ ^' q ^ dy. 
Jr* Jr* 

The phase function is given by 

©(», q,p) = -p-(y-q) + ||y - q\ 2 + s {y)- 

Clearly, 3m > and 3m = if and only if y = q. The derivatives of with 
respect to y arc 

d y G = p + d y S a (y) + i{y - q), 

d 2 y e = dls (y) + iL 

Hence, the first derivative vanishes only when y = q and p — d y So(y), and 
dct<9^0 7^ 0. Therefore, we can apply stationary phase approximation with com- 
plex phase (see for example [7]) to conclude, for (q,p) £ K 2d , 

/ c{q)\p\e-c p - {v - q) -^ y - q] > 2 ' A ± (y)ei So{y) dy 
Jm d 

= f c(y)\d y S (y)\e-^< y - q '>-^ y - q fA ± (y)ei Sa Mdy + O(e). 

JR d 
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Therefore, 
(3.23) 

2 d/ 

d t u(0,x) -- 



(2ire) 3d / 2 
2<i/2 



[ -c(y)\d y S (y)\ei^ - x 'y^u +fi (y)dydpdq 
[ -c{y)\d y S {y)\e^-^^u_ fi {y)dydpdq + O{l). 



(27T£) 3d / 2 

Substitute (|3.22l) into (I3.23|) and use Lemma T3.ll then 

d t u(0, x) = -B (x)ei s °W, 
e 

which agrees with (|1.2[) . 

3.2. Derivation of the evolution equation of a±. In order to derive the evo- 
lution equation for the weight function a, we carry out the asymptotic analysis of 
the wave equation (jl.ip using the ansatz p. II) in this section. As the equation (jl.ip 
is linear, we can deal with the two branches separately. In the following, we only 
deal with the "+" branch that corresponds to H + , and the other is completely 
analogous. For simplicity, we drop the subscript "+" in the notations. 

Substituting (|3.1j) into the equation (|1.1|) (keeping only the "+" branch) gives 

d 2 t u = {2 J )3d/2 J R3d ( d ? a + 2 ~ £ d t ad ^ + ~ £ ad t^ ~ ^a(d t <S>) 2 )e^^u dydpdq, 
and 

Am = \ 3d/2 I \ {d x <b ■ d x $>))ae^ It Ua dy dp dq. 

Squaring both sides of (|3.15[) yields 

(3.24) (d t <S>) 2 = c 2 \P\ 2 + ((x - Q) ■ (\P\d Q c + )) 2 

+ 2c\P\(x-Q)-(\P\d Q c+zc^ 
Differentiating (|3.15l) with respect to t, one has 
d 2 f> = -d t (\P\c) + d t Q ■ (\P\d Q c 

(3.25) - (x - Q) ■ (d Q c^j^- + d 2 Q c ■ d t Q\P\ 



P 

)C + l-rr=rrC 



d t p d p -d t p p 

»C-r^r - icP h i-t^Oqc ■ a t Q 



\p\ p ■ 

We simplify the last equation using (12.41) . 

d 2 § = cP ■ d Q c + ic 2 



^3 26) -(x — Q)- \-d Q cP ■ d Q c + cd Q c ■ P 

,Pd Q c 



icdnc + 2icP- 
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Taking derivatives with respect to x produces 

(3.27) d x <S> = P + i(x- Q), 

(3.28) c* x $ • 3 X $ = \P\ 2 + 2iP .(x-Q)-\x- Q\ 2 , 
and 

(3.29) A$ = di. 
We next expand c{x) around the point Q, 

(3.30) c{x) =c + d Q c-(x-Q) + -(x - Q) ■ d 2 Q c{x - Q) + 0(\x - Q\) 3 , 
and 

(3.31) c 2 (x) = c 2 + 2cd Q c -(x-Q) + (d Q c ■ (x - Q)) 2 

+ c(x - Q) ■ d 2 Q c{x -Q) + 0(\x - Q\) 3 . 

The terms c, 8qc and DqC on the right hand sides are all evaluated at Q. 

Substituting all the above into the wave equation (|1.1|) and keeping only the 
leading order terms give 

2-d t a(-c\P\)u + -a(cP ■ d Q c + ic 2 )u 

e e 

- ^a(2c(x - Q) ■ (\P\ 2 d Q c + icP) + ((x - Q) ■ (\P\d Q c + icP /\P\)) 2 )u 

- c 2 (-^-ad - ^aP -(x-Q) + -^a\x - Q\ 2 ^ju 
+ ^acdqc ■ (x - Q){\P\ 2 + 2iP ■ (x - Q))u 

+ ^a\P\ 2 ((d Q c ■ (x - Q)) 2 u + c(x - Q) ■ d 2 Q c{x - Q))u ~ 0(1). 
After reorganizing the terms, we get 

(3.32) 2-c\P\d t au ~ -a(cP ■ d Q c -(d - l)c 2 i)u - \a(x - Q) ■ Mix - Q)u, 
where 

(3.33) M = (\P\8qc - icP/\P\) ® (\P\8qc - icP/\P\) + c 2 I 

- \P\ 2 d Q c®d Q c-\P\ 2 cd 2 Q c. 

Lemma l3~2l shows that 

(3.34) a(x - Q) ■ M(x - Q)u ~ eatr(Z" 1 d z QM)u + 0(e 2 ). 
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Therefore, to the leading order, we obtain the evolution equation of a, 
(3.35) ft.-K*.,*-^ 

+ |tr(z-'S,Q(2^«,a c- ^(^& - I) - AP\S' QC 

Notice that 

dZ fdQ dP\ ( P 

and 



PI 

By using the fact that (|3.34p has a quadratic form, one has 

tr (V^Q^ ® 8 Q ^j = tr f Z-^.Q^ ® P 

Hence (13.35)) can be reformulated as 

da P a / i dZ 

— = a- — r • <9qc + - tr Z' 1 

dt \P\ Q 2 \ dt 

This completes the asymptotic derivation. We remark that in the case of time 
dependent Schrodingcr equation, the asymptotics have been made rigorous in [25, 

Eg. 

4. Numerical examples 

In this section, we give both one and two dimensional numerical examples to 
justify the accuracy of the frozen Gaussian approximation (FGA). Without loss of 
generality, we only consider the wave propagation determined by the "+" branch 
of (|2~Tj) which implies that B (x) = -ic(x)\W x Sq(x)\A (x) in (fl~2]) . 

4.1. One dimension. Using one-dimensional examples in this section, we compare 
FGA with the Gaussian beam method (GBM) in both the accuracy and the per- 
formance when beams spread in GBM. We denote the solution of GBM as u GBM , 
and summarize its discrete numerical formulation (only the "+" branch) as follows 
for readers' convenience ([14 U'241I29] ). 



3=1 

x exp i 



i GBM (t,x) = £ (^) 2 r e (\x - yj \)A(t, yj ) 

(i(S(t, Vj ) + £(t, yj )(x - y 3 ) + M(t, y 3 )(x - Vj f /2)) Sy , 
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and yj, £, S, M, A satisfy 



-^=c{y j )-^-, Vj(0)=V } , 

dt 
dS 
~dt 



— = 0, S(0)=S ( yj ), 



^ = -2d yj c(y 3 ) j|m - d 2 y] c( yj )\t\, M (0) = d 2 y . S (yj )+t, 
dt 2 7 ' !h> < 



where rg is the cutoff function, j/q's are the equidistant mesh points, Sy is the mesh 
size and N yo is the total number of the beams initially centered at y J Q . 

Example 4.1. The wave speed is c(x) = x 2 . The initial conditions are 
uo = cxp(-100(a; - 0.5) 2 ) exp (y) , 

Q 

d t u = — — exp(-100(x - 0.5) 2 ) exp (-j . 



The final time is T = 0.5. We plot the real part of the wave field obtained by 
FGA compared with the true solution in Figure [2] for e = 1/64, 1/128, 1/256. The 
true solution is computed by the finite difference method using the mesh size of 
Sx = 1/2 12 and the time step of St = 1/2 18 for domain [0,2]. Table [T] shows the 
i°° and I 2 errors of both the FGA solution u FGA and the GBM solution u GBM . 
The convergence orders in e of i°° and t 2 norms are 1.08 and 1.17 separately for 
FGA, and 0.54 and 0.57 for GBM. We observe a better accuracy order of FGA than 
GBM. 

We choose St = 1/2 11 for solving the ODEs and Sx = 1/2 12 to construct the 
final solution in both FGA and GBM. In FGA, we take Sq = 5p = Sy = 1/2 7 , N q = 
128, N p = 45 for e = 1/128, 1/256 and Sq = Sp = Sy = 1/2 5 , N q = 32, N p = 33 
for e = 1/64. In GBM, we take Sy = 1/2 7 , N ya = 128 for e = 1/128, 1/256 and 
Sy = 1/2 5 , N yo = 32 for e = 1/64. 

We remark that in this example the mesh sizes of p and q have been taken 
very small and N p large enough to make sure that the error of FGA mostly comes 
from asymptotic expansion, but not from initial value decomposition, numerical 
integration of ODEs and so on. Such a choice of fine mesh is not necessary for the 
accuracy of FGA, as one can see in Example 14.31 
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£ 


1/2 6 




1/2 7 




1/2 8 




lltt- 


u FGA \\ e oo 


1.12 x 10" 


-l 


6.18 x 10- 


-2 


2.51 x 10" 


-2 


ll«- 


u FGA \\ P 


6.05 x 10- 


-2 


2.96 x 10" 


-2 


1.19 x 10- 


-2 


||u- 


u GBM |U- 


7.15 x 10" 


-1 


5.08 x 10" 


-1 


3.36 x 10" 


-1 


ll«- 


U GBM \\ e 2 


3.26 x 10" 


-1 


2.28 x 10" 


-1 


1.47 x 10" 


-1 



Example 4.2. The wave speed is c(x) 

(a; -0.55) 



Uo = exp l v - ~ 2g - 

ix 1 ( (x — 0.55) 2 \ (ix\ 
9 tU0 = -_exp( )exp(-j 



We use this example to illustrate the performances of FGA and GBM when the 
beams spread in GBM. The final time is T — 1.0 and e = 1/256. Remark that 
the initial condition is chosen as a single beam on purpose so that one can apply 
GBM without introducing any initial errors. The true solution is provided by the 
finite difference method using Sx = 1/2 11 and St = 1/2 17 for domain [0, 4]. We take 
Sq = 5p = Sy = 1/2 7 , N q = 128, N p = 45 in FGA to make sure that the error in 
the initial value decomposition of FGA is very small. The time step is St — 1/2 10 
for solving the ODEs and the mesh size is chosen as Sx = 1/2 11 to construct the 
final solution in both FGA and GBM. 

Figure [3] compares the amplitudes of the wave field computed by FGA and GBM, 
and the true solution. One can see that the beam has spread severely in GBM. The 
results confirm that FGA has a good performance even when the beam spreads, 
while GBM does not. Moreover, it does not help improving the accuracy if one 
uses more Gaussian beams to approximate the initial condition in GBM as shown 
in Figure HI where N yo = 128 beams are used initially and Syo = 1/2 7 . Remark 
that GBM can still give good approximation around beam center where Taylor 
expansion does not introduce large errors. This can be seen around x = 1.2 in 
Figure [3] 

4.2. Two dimension. 

Example 4.3. The wave speed is c(x\,X2) = 1. The initial conditions are 
uq = exp(— lQQ(xi + a;|)) exp x i + cos(2a;2))^ , 



d t u = --y 1 + 4sin 2 (2iE2) exp(-100(a; 2 + x\)) exp y-{-xi + cos(2x 2 )) 

This example presents the cusp caustics shown in Figure [5] The final time 
is T = 1.0. The true solution is given by the spectral method using the mesh 
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(b)e 




(c)e 
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E 0.04 

Z 0.02 


0.2 

1 

" 128 
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2 0.08 

g 0.06 
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E 0.04 
Z 0.02 




0.4 



0.6 O.f 
X 



0.2 



0.4 



0.6 O.f 
X 



1 

256 



Figure 2. Example 14 .!( the comparison of the true solution (solid 
line) and the solution by FGA (dashed line). Left: the real part of 
wave field; right: the errors between them. 



8x\ = 8x 2 = 1/512 for domain [—1.5,0.5] x [—1,1]. We take dqi = 8q 2 = Spx — 
S'P2 = Syi = Sy 2 = 1/32, N q = 32, N p = 8 in FGA, and use Sx x = Sx 2 = 1/128 to 
reconstruct the solution. Figure [6] compares the wave amplitude of the true solution 
and the one by FGA for e = 1/128 and 1/256. The £°° and £ 2 errors of the wave 
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Figure 3. Example l4.21 the comparison of the true solution (solid 
line) , the solution by FGA (dashed line) and the solution by GBM 
(dots) for e = ^g. Left: the amplitude of wave field; right: the 
error between them (dashed line for FGA, dots for GBM). 




0.5 1 1.5 2 2.5 3 3.5 
X 



0.8 




0.5 1 1.5 2 2.5 3 3.5 
X 



Figure 4. Example l4.2i the comparison of the true solution (solid 
line) and the solution by GBM using multiple Gaussian initial rep- 
resentation (dots) for e = ^|g. Left: the amplitude of wave field; 
right: the error between them. 



amplitude are 1.98 x 1CT 1 and 4.42 x 1(T 2 for e = 1/128, and 1.07 x 1CT 1 and 
2.20 x 10~ 2 for e = 1/256. This shows a linear convergence in e of the method. 



5. Discussion and Conclusion 

We first briefly compare the efficiency of frozen Gaussian approximation (FGA) 
with the Gaussian beam method (GBM). GBM uses only one Gaussian function for 
each grid point in physical space, while FGA requires more Gaussians per grid point 
with different initial momentum to capture the behavior of focusing or spreading 
of the solution. However, the stationary phase approximation suggests that the 
number of Gaussians is only increased by a small constant multiple of the number 
of those used in GBM. In addition, in GBM one has to solve the Riccati equation, 
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the cusp caustic. 

which is a coupled nonlinear ODE system in high dimension, to get the dynamics 
of the Hessian matrix for each Gaussian, while in FGA the Hessian matrix is de- 
termined initially and has no dynamics. Therefore, the overall efficiency of FGA is 
comparable to GBM. 

Admittedly, higher order GBM gives better asymptotic accuracy, and only re- 
quires solving a constant number of additional ODEs as in FGA. The numerical 
cost of higher order GBM is comparable to FGA. However, higher order GBM has 
its drawbacks: The imaginary part of higher order (larger than two) tensor func- 
tion dose not preserve positive definiteness in time evolution, which may destroy 
the decay property of the ansatz of higher order GBM . This is even more severe 
when beams spread. Moreover, the ODEs in higher order GBM arc in the form of 
coupled nonlinear system in high dimension. It raises numerical difficulty caused by 
stability issues. We also note that, the numerical integration of ODEs in FGA can 
be easily parallelized since the Hamiltonian flow (|2.4j) is independent for different 
initial (q,p), while it is not so trivial for higher order tensors in GBM. 

From the accuracy point of view, our numerical examples show that first or- 
der FGA method has asymptotic accuracy 0(e). The existing rigorous analysis 
([2[2j[T6]) proves that the fc-th order GBM has an accuracy of 0(e k ^ 2 ). Hence, 
at the first order, FGA has better asymptotic accuracy than GBM. We note that, 
however, there has been numerical evidence presenting 0(e) asymptotic accuracy 
order for first order GBM, for example in [9,11,16, 18 . This phenomenon is usually 
attributed to error cancellation between different beams. To the best of our knowl- 
edge, the mechanism of error cancellation in GBM has not been systematically 
understood yet. 

With the the gain of halfth order in asymptotic accuracy due to cancellation, 
the first order GBM has the same accuracy order as FGA (of course GBM still 
loses accuracy when beams spread). Remark that the gain in asymptotic accuracy 
order depends on the choice of norm. For example, the first order GBM has a 
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(c) Errors 

Figure 6. Example I4.31 the comparison of the true solution and 
the solution by FGA. Left: wave amplitude of e — \ right: wave 
amplitude of e = • 



halfth order convergence in £°° norm, first order convergence in £ 2 norm and 3/2-th 
order convergence in i 1 norm in Example 1 of 9 . Moreover, the error cancellation 
seems not to be easily observed in numerics unless e is very small. For instance, the 
convergence order of GBM in Example 14. II is only a bit better than 1/2 for s up to 
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1/256. While in FGA, we numerically observe the first order asymptotic accuracy 
in both £ 2 and t°° norms. 

Actually the accuracy of FGA can also be understood from a viewpoint of error 
cancellation. Note that the equalities (|3.9p . (|3.10l) and (|3.11l) in Lemma 13721 play the 
role of determining the accuracy of FGA. In (|3.9[) . the term x — Q is of order 0(y/e) 
due to the Gaussian factor, but after integration with respect to q and p, which is 
similar to the beam summation in GBM, it becomes 0(e). Similar improvement 
of order also happens in (j3. 10[) and p. lip . Integration by parts along with (|3.14p 
explains the mechanism of this type of error cancellation. 

We conclude the paper as follows. In this work, we propose the frozen Gaussian 
approximation (FGA) for computation of high frequency wave propagation, mo- 
tivated by the Herman-Kluk propagator in chemistry literature. This method is 
based on asymptotic analysis and constructs the solution using Gaussian functions 
with fixed widths that live on the phase plane. It not only provides an accurate 
asymptotic solution in the presence of caustics, but also resolves the problem in 
the Gaussian beam method (GBM) when beams spread. These merits are justified 
by numerical examples. Additionally, numerical examples also show that FGA ex- 
hibits better asymptotic accuracy than GBM. These advantages make FGA quite 
competitive for computing high frequency wave propagation. 

For the purpose of presenting the idea simply and clearly, we only describe the 
method for the linear scalar wave equation using leading order approximation. The 
method can be generalized for solving other hyperbolic equations and systems with 
a character of high frequency. The higher order approximation can also be de- 
rived. Since the method is of Lagrangian type, the issue of divergence still remains, 
which will be resolved in an Eulerian framework. We present these results in the 
subsequent paper |17) . 
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